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Abstract 

A compact difference scheme derived in [3] for treating the equilibrium 
equations of elasticity is studied. The scheme turns out to be inconsistent 
and unstable. A multigrid method which takes into account these properties is 
described. The solution of the discrete equations, up to the level of 
discretization errors, is obtained by this method in just two multigrid 
cycles. 
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1 . INTRODUCTION 


In this paper, we study a compact finite difference scheme for the 
equilibrium equations of elasticity derived in [3]. We focus here on the two- 
dimensional case only. 

We begin in Section 2 by the derivation of the scheme for a general 
source terra in the elasticity equations (i.e., a non-equilibrium case). This 
is needed later when a multigrid method is considered. 

In the first step, equations for displacements and stresses in a cell are 
obtained. They consist of equations which approximate the equilibrium 
equations and those which represent some single-valuedness of displacements in 
cell centers. The second step is a process of elimination of stresses. This 
results in a set of equations involving displacements only. The equations 
divide the set of grid points into two dispoint sets; on each of them 
different equations are given. 

Section 3 deals with the inconsistency of the resulting scheme. A Taylor 
expansion of the different terms for either of the sets of equations shows 
that both are inconsistent with the equations of elasticity. However, a 
closer look reveals that consistency-in-the-average exists. That is, the sum 
of the equations in a cell shows the desired consistency. This fact is used 
later when a multigrid method for that scheme is derived. 

Section 4 deals with the instability of the scheme. It is shown by means 
of Fourier analysis that the interior equations admit highly oscillatory 
solutions for the homogeneous problem. This means that there exist non-smooth 
components that will not affect the residuals almost at all. This is a 
numerical instability since it does not have a differential analogue. Schemes 
with such a property have been studied in [2] and are referred to as quasi- 
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elllptic schemes. In particular, a multigrid method for such schemes is 
described there. 

Sections 5 and 6 describe the multigrid ideas and their implementation in 
the present context. A standard Gauss-Seidel relaxation was used for relaxing 
interior equations, while a modification is introduced when traction boundary 
conditions are relaxed. The coarsening used was done in a compatible way to 
the discretization. An averaging operator is applied to grid functions before 
interpolating them to finer levels. This is done in order to remove unstable 
components in the approximation. The resulting full multigrid algorithm 
(nested iteration) solves the discretized equation, up to the level of 
discretization errors in just two cycles. Thus a very efficient method of 
solution is obtained. 


2. COMPACT SCHEME 

In two dimensions, the stresses t = (Tj, * 2 ) are given in terras of 
the displacements u = (uj, U 2 ) as 


T 11 = ?3 Xl u l + n9 x 2 u 2 


(2.1a) 


t 22 = ^“l + C3 x 2 u 2 


(2.1b) 


t 21 ' t 12 ‘ 2 0(S * 2 U 1 + 3 Xj U 2>' 


(2.1c) 


The parameters £ , ti , a are given in terms of Young's modulus E and 


Poisson's ratio v by 



? * 

( 1 -/) 


h = 


Ev 

? * 

<1-0 


E 

Tl+vJ * 


The equations of elasticity in terms of the stresses are 


3 t. . + 3 x._ - f . 

x^ 11 12 1 


(2.2a) 


^^21 * 3 x 2 T 22 f 2* 


(2.2b) 


where in equilibrium f^ = ? 2 = 0. We assume here a nonequilibrium case, 
i.e., f^ ? 2 are not zero. This is needed later when coarse grid equations 
are considered. 

A compact scheme for square cells has been derived in [3]. It is assumed 
that displacements and stresses are given on cell midfaces as shown in Figure 
1 . The scheme is given by 


\ T n ' + 


(2.3a) 


“ X2 t 22 ' n5 Xl u l + ?S * 2 u 2 


(2.3b) 


IJ Xj T 21 ‘ "x 2 t 12 " 2 + 5 x 1 u 2 ) 


(2.3c) 


6 T. . + 6 f 10 = f , 

x^ 11 x 2 12 1 


(2.3d) 


6 x 1 T 21 + 6 x 2 T 22 f 2 


(2.3e) 


2 2 
Vi u, - kh 5 t, t = y u. - kh 6 t, 0 
H x^ 1 x^ 11 x^ 1 x^ 12 


(2.3f ) 
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y Xi u 2 - kh 2 5 Xi t 21 = ^u 2 - kh 2 6 x2 t 22 . 


(2.3g) 


where fj, f£ are midcell values of fj, f 2 , respectively, and the operators 


y , 6 (i = 1, 2) are defined by 
X i x i 


U x ^ x ^ = ^ X 1 + h ’ x 2 ^ + ^( x i“ h > x 2 ^ / 2 


5 x = W x i + h l x 2^ " ^ x i~ 


and similarly y , 6 . The parameter k is arbitrary. 

2 X 2 


U, T 



u, T 


Figure 1 

In actual computation the stresses are eliminated from the equations and 
one gets equations involving the displacement only. We will redo here the 
elimination of the stresses, for the case equation (2.3d) through (2.3e) has 
nonzero right-hand sides. 

Using (2.3d), (2.3e) in (2.3f), and (2.3g) respectively, we get 


h6 t 
x i 


ii 


2kh ^ u x. 


\ )u l + 2 f l 


(2.4a) 
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hs x 2 T 22 ■ TOT - “x 2 >u 2 + T V 


(2.4b) 


Also we have from (2.4) and (2.3d), (2.3e) the relations 


hS x, T 21 ’ 2® ( “x 2 - \ )a 2 + l f 2 


(2.4c) 


h6 x 2 T 12 2kh (y x 2 “ y Xl )u l + 2 f l* 


(2.4d) 




Figure 2 


The idea of elimination is to set equal the values of t that are computed 
from two neighboring cells. That is, 


[«** " y x " K + h6 x >] 

X 1 1 X 1 X 1 


T 1 = 0 


(2.6a) 




-6- 


g 

[“x 2 • “x 2 ■ <hS x 2 + h?S 2 ) l t 2 ’ 0 

where y R , , w® (i - 1,2) denote averages 

i i i i 

and bottom cells, respectively, and correspondingly h6 x 
denote differences (divided by 2) on right, left, top, and bottom cells 
respectively (see Figure 2). 

Upon inserting (2.3a-c) and (2.4) into (2.6), we get the following 
equations: 


(2.6b) 

on right, left, top, 
, h5 L , h5 T , hS B 

7 V 7 V 7 V 


S(6 R - 6 L )u. + n(<$ R 

x l x j 1 x 2 


S x 2 >“2 - 


1 / R R % , / L L \ 

2kh ( \ ' %> + ( \- \> 

I Z • 1 z 


'] , 


I (f R + fj) = o 


(2.7a) 


1 y.R pL V 1 f 

2 a<S x 2 ‘ 5 x 2 )u 1 + 2 0<{ Xj 


^>“2 ' 


1 ( R 
2kh ''V 


R \ / L L v 

u ) + (u - u ) 


x. 


Xo 


u 2 * 
(2.7b) 


| <£* + f 2 > - 0 


I a({ x 2 " S x 2 )u l + 7 " 8 x , )u 2 - 


1 ( T T * , , B B v 

2kh ( "x 2 ‘ + <u x 2 ‘ %> 


u, - 


h 

2 



+ fj) - 0 


(2.7c) 



(2.7d) 


"<**, - 5 ^"I + ' W x 2 


S “ 2 )u 2 ' 


1 , T T B B 

Ikh ‘ + (l, x 1 - “x 2 > 


Uo “ 


| ( f 2 + f 2> 


0. 


Note that equations (2,7a) and (2.7b) are given on P-points, while (2.7c) 
and (2.7d) are given on Q-points (see Figure 3). 



Figure 3 

Since it is natural to expect that some of the boundary conditions will 
be given in terms of stresses, we have to express the stresses at the boundary 
in terms of the displacement. That is, 


t 22 B ■ "‘xj”! + 5{ X„ U 2 * TO (u x. - u 2 * 1 f 2 


J Un 

x 2 2 


(2.8a) 
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T 12 2 ° (6 x 2 U 1 + 6 Xj U 2 J * 2kh (y x 2 “ “x^ U 1 * 2 f l 


(2.8b) 


T 12 L = ?5 Xl U l + n6 x 2 U 2 * 20T ” y x 2 ) u l ± I f l 


(2.8c) 


T n L = 7 0(6 x 2 U 1 + 6 x 1 U 2 ) * 2^ (y x 2 " y Xl > “2 * I £ J 


(2.8d) 


In the next sections we assume the parameter k to have the value — . 


3. INCONSISTENCY AND IMPLICATIONS 

The equations of elasticity, in terms of displacements, are, 


2 2 12 1 ~ 
£3 u. + ti3 u 0 + 7 T a 3 u. + ^ od u n = f . 
* 1 x 2 2 2 x 2 1 2 x i x 2 2 1 


4- cr3 u, + i a3 2 u 0 + n3 u, + ?3 2 u 0 = f 0 , 
2 x i x 2 1 2 2 x l x 2 1 x 2 2 2 


(3.1a) 

(3.1b) 


The compact scheme (2.7) is not consistent with this equation in the usual 
sense. In fact, if we expand the different terms in (2.7) in a Taylor 
expansion, we find the following: The equations at P-points are consistent 

with 


i Au l + - | )3 X, U 1 + ^X.x "2 = 7 ? 1 


1 


k l~2 


a ~ . a . 1 -x 

7 S Xl x 2 U l + T 4u 2 ’ 2 f 2- 


(3.2a) 

(3.2b) 


and the Q-equations are consistent with 



(3.3a) 


a . , 1 . 1 y 

T 4u l + 2 “^x^ ■ T f 1 


+ T 4u 2 + (5 ‘ f> 8 x 2 u 2 ■ 1 V 


(3.3b) 


Neither (3.2) nor (3.3) are the equilibrium equation of elasticity. However, 
by summing (3.2) with (3.2) and using the definition of £, ti, a, we get back 
(3.1). That is, we have consistency in an averaged sense. 

The fact that we have only consistency in-the-average is important when 
coarsening has taken place in a multigrid process. Residuals that are 
transferred to coarser levels should better be related to equations which are 
consistent with the differential equation. 

This can be achieved for equations (2.7) if we sum the proper equations 
to ensure consistency with (3.1). That is, by summing (2.7a) with (2.7d) and 
(2.7b) with (2.7c), consistency with (3.1a) and (3.1b), respectively, is 
obtained. 

If boundary conditions are given in terras of stresses, residual transfer 
should be done carefully. By looking at (2.8) we see that the discretization 
is such that the boundary condition has contribution from the right hand side 
of the interior equation. Residual transfer to coarser levels should maintain 
this. 


4. INSTABILITY 

Another important property of the scheme (2.7) is its instability . That 
is, there are highly oscillatory displacements which satisfy the interior 
equations with zero right hand side. This means that large changes in these 
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components, or components close to them, will not affect the residual almost 
at all. Hence, a small change in the equations can introduce large changes in 
the solution in these unstable components. This is a numerical instability, 
since a correspondingly large change in the differential solution cannot 
occur. 

We show below how to find the unstable components by means of Fourier 
analysis. We consider the homogeneous equation and write it in the form 



where U p (Uq) denotes the uj values on P(Q) points and similarly V p 
(Vq) denotes U 2 values on P(Q) points. Consider displacements of the 
form 



where 0 - (0. ^ , 8 2 ) > H = 2h (the size of a cell), and 1 0 | < tt, 

1 0 | = max(|0jj, 1 0 2 | ) • For this choice of displacements we have 


where 
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where (0) is a 4x4 matrix of functions depending on JK By 

looking at (0) we can examine the stability of our scheme. If there 

exists a 0_ * 0 such that det (0) = 0, it means that there are high 
frequency components which solve the homogeneous equation. For our scheme 

A 

det <£ (O) = 0 


implies 0_ = (0,0) or The unstable Fourier component is 

therefore 0 = (tt, tt). Its amplitude is obtained by solving the equation 



for a nontrivial solution. This gives us 
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This corresponds to the following displacements (in a cell). 



The computation above shows that these are the only unstable 
components. Hence, this scheme satisfies the definition of quasi-ellipticity 
given [2]. 

Since det (0) = 0 for some _0 * 0, in an infinite space, or 

under periodic boundary conditions, there exists a highly oscillatory function 
h 

v (x) = A exp (i0_*x/H) which satisfies the homogeneous equation 
h h 

L v n = 0, Hence, the solution, unlike the corresponding differential 

solution, is not unique (up to an additive constant); it contains an 

undetermined highly oscillatory component. Similarly in any bounded domain 

h h 

with any boundary conditions, functions W n (x) close to v n (x) 
h h 

(e.g,, W n = <|>i v n + » <t>j being smooth) exist which satisfy the boundary 

L V -L 

condition and for which L n W n is everywhere small. Such W , therefore, forms 
an unstable mode: a small change in the equation can introduce a large change 

proportional to W • 

This is a kind of numerical instability, since a correspondingly large 
change in the differential solution cannot occur. The numerical instability 
need not hurt much: If the differential system is LU = f and the discrete 

system is L^U* 1 = f* 1 , all one has to do is to define F* 1 = I^F, say, through 
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an averaging operator which liquidates the unstable modes. Another way to 
remove the instability is by averaging the solution, that is, by replacing the 
computed solution U* 1 by S^U* 1 , where is an averaging operator which 

removes all the unstable components but retains the accuracy of smooth 

components. 

When derivatives are calculated, much greater loss of accuracy can occur 
than in the solution itself. The averaging process discussed above reduces 
this inaccuracy in derivatives. 

In the problem we treat here, the stresses, as computed from the 

displacements, are the important physical quantities. Since they involve 
derivative of displacements, we might expect to see degradation of accuracy as 
a result of the instability. However, the unstable components for the scheme 
discussed here are such that they produce zero stresses. That is, no loss of 
accuracy in the stresses occurs as a result of the instability. 

The mentioned instability has strong implications on the multigrid method 
to be used. Usual multigrid solvers yield poor asymptotic rates when applied 
to quasi-elliptic schemes. The reason is simple: slow to converge are the 

unstable modes. They cannot converge by coarse-grid correction, since they 
are high-frequency modes, essentially invisible on coarser levels. Neither 
can they significantly converge by any type of local relaxation since these 

unstable modes show a very small residual function (compared with residuals 

shown by other modes with comparable amplitude) and the correction introduced 
by relaxation is proportional to the size of the residuals (see [1]). The 
smoothing factor for such schemes is 1, and it is achieved at the unstable 


modes • 
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The poor asymptotic convergence is not important. The modes which are 
slow to converge are exactly those unstable ones for which algebraic 
convergence is not really desired, their amplitudes in the algebraic solution 
being unrelated to their amplitude in the differential equation. The only 
concern is that these amplitudes will remain suitably small. 


5. PRACTICAL IMPLEMENTATION 

In this section we describe the elements of a multigrid procedure defined 
in Section 6. 

Relaxation 

As a relaxation we have used Gauss-Seidel for the interior equation, 

where the P-points are relaxed first followed by the Q-points. 

Relaxing the traction boundary conditions is done slightly differently 

(see [1], Section 5.3). Instead of performing Gauss-Seidel for the traction 

boundary condition BU = g, we perform Gauss-Seidel on the equation 
2 2 

3 ^ d 3 

— _ BU = — =• g, where ■=— is derivative tangential to the boundary. 
ds Z 3s Z dS 

Practically, it means that instead of satisfying a given equation at a 
boundary point, we only change it such that its error is the average of the 
errors at neighboring boundary points. This relaxation preserves the 
smoothness of the interior approximation, unlike the straightforward Gauss- 

y 


Seidel 
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Coarsening 

The coarsening method we have used here is referred in [1] as compatible 
coarsening. That is, the coarsening procedure is analogous to the fine-grid 
discretization. 



O Coarse grid points 
• Fine grid points 


Figure 5 


Each 2x2 fine-grid cell forms a coarse grid cell [see figure 5]. Coarse 
grid equations are defined using formulas (2.7) where f^, are replaced by 
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residuals of finer levels, which are consistent with differential residuals 
(as mentioned in section 4). 

Interpolation 

In interpolating corrections from a coarse to a fine grid, we used the 
following procedure: 

(a) Define values at midcells and at vertices on coarse level by linear 
interpolation. 

(b) Linearaly interpolate corrections to the fine grid, using midcells 
and vertex values defined in (a), and original coarse-grid function. 

Averaging 

An averaging operator was applied to corrections before interpolating 
them. This was done in order not to introduce high frequency errors which are 
unstable. The averaging was such that it killed the unstable mode but 
retained the accuracy of the scheme. It is given by the following formula: 

(sV)(x) = I [4u^(x,y) + + h, y + h) + 

u ± (x + hj , y - h) + u ± (x - h, y + h) + u^x - h, y - h)]. 


(i = 1,2). 
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6. FMG SOLUTION TO TRUNCATION LEVEL 

Since the raultigrid cycling is inefficient in reducing unstable mode 
errors, the multigrid solver should take care not to start with an initial 
solution which contains large amplitudes of such errors. The overall initial 
error in unstable modes should be smaller than the overall truncation error. 
This is easily obtained by taking a first approximation from a coarser grid, 
employing interpolation of suitable order. The usual n Full multigrid" (FMG, 
also called "nested iteration") algorithm can therefore be used, with slight 
modifications described in the following. For a flowchart, and a 
detailed discussion of FMG algorithms and the order of the first 

interpolation, see Secs. 1.6 and 7 in [1]. We describe here the Correction 
Scheme (CS) version of the algorithm since our problem is linear, and issues 
of local refinement are not discussed here. 

6.1 Multigrid cycle. 

Suppose a sequence of grids is given with meshsizes h^ (k = 1, 2, 

3,...) where h^ + ^ = h^/2. On the h^ grid the discrete equations have the 
form 

L k U k « F k (6.1) 

where L k approximates L k+ ^. Given u k , an approximate solution to 

(6.1), the multigrid cycle MG for producing an improved approximation u k 


u k +MG(k, u k , F k ) 


( 6 . 2 ) 
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is recursively defined as follows: 

If k = 1 solve (6.1) by any direct or iterative method, yielding the 
k 

final result u^. Otherwise do (A) through (D). 


(A) Perform 
approximation 



relaxation sweeps on (6.1), 


resulting in a new 


k— 1 

(B) Starting with Uq =0, make y successive cycles 

k™ 1 f i i k- 1 - k— 1 1 _ k T k-k ^ ^ / • i \ 

Uj -«-MG(k - 1, Uj j, I fc (F - L u )), (j = 1,..., y) 

k— 1 

where 1^ is a transfer ("reduction”) of residuals from grid 
\ to grid h k-1 . 

(C) Calculate u^ = u^ + 1^ , * 

k-1 y * 

interpolation ("prolongation") from grid 
k-1 

S is a suitable averaging operator. 

~k 

(D) Perform V 2 relaxation sweeps on (6.1), starting with u and 
yielding the final result u^ • 

The cycle with y = 1 is called V cycle or V( ), and the one 

y = 2 is called W cycle or W( Vp ). 


k 

where I, , is a suitable 
k-1 

h k-l t0 grid h k» and 


with 
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6.2 Full Multigrid (FMG). 

The N-FMG is an algorithm for calculating an approximate solution 

u k = FMG(k, F k , N) (6.4) 


to equation (6.1), defined recursively as -the following two successive steps. 

k 

(a) Calculating a first approximation u^: 

v 

if k = 1, put u q = 0. Otherwise, put 

u k = n k _ 1 FMG(k - 1, l£ -1 F k , N), (6.5) 

k 

where is an interpolation operator from grid to 

k-1 

grid h^, and 1^ is a transfer from grid k to grid k-1. 

k 

In our experiments, the interpolation operator II was the same 

as the one for interpolating corrections (defined in Section 5). 


(b) Improve the first approximation by N successive MG cycles 
u k +MG(k, u k _ x , F k ), (j = 1,..., N) 


as defined in Sec. 6.1 
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7. RESULTS AND DISCUSSION 

A domain {(x, y):0<x<l, 0<y<l} is considered. Boundary 
conditions are given in terras of displacements on the two boundaries x = 0, 
and x « 1. Stresses are given on the boundaries y = 0, y = 1. Source terras 
are used such that the exact solution is given by 

u l = u 2 ~ sin (»2(x + 2y - 2)). 

The finest level uses a mesh size of 1/32, and 4 levels were used in the 
raultigrid process. In the notation of Section 6 the following parameters were 
used: y ■ 2, = 3, V 2 - 3, N - 4, The table below shows the dynamic L 2 - 

norm of the residuals on the currently finest level, as well as the error in 
the approximation at that stage. It is seen that after two cycles the problem 
is solved to the level of discretization errors. We show results for N - 4 
not because it is really needed, but in order to show that indeed, N = 2 is 
enough to obtain approximate solution whose errors are below the level of 
discretization errors. Observe that the asymptotic convergence rate, as 
predicted, is not very fast. It is related to the unstable components. This 
should not bother us since we do not want the unstable component to 
converge. All we need from these components is to have errors below the level 
of discretization errors. 

The results clearly show that a very efficient method for dealing with 
compact scheme is obtained in spite of the fact that such schemes have 
instability and inconsistency properties. 



- 21 - 


Level # 

Cycle // 

IIResiduals^ 

»u h - vn 

2 

5 

.505(-4) 

•406(-3) 

3 

1 

.482(-2) 

.429(-4) 


2 

.109(-2) 

•274(-4) 


3 

•279(-3) 

.252(-4) 


4 

•861 (—4) 

.247(-4) 

4 

1 

• 132(— 2) 

• 67 6(— 5) 


2 

•520(-3) 

•558(-5) 


3 

•288(— 3) 

•558(-5) 


4 

• 1 62(— 3 ) 

.558(-5) 
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